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The propagation of a wave-packet in a nonlinear disordered medium exhibits interesting dynam- 
ics. Here, we present an analysis based on the nonlinear Schrodinger equation (Gross-Pitaevskii 
equation) . This problem is directly connected to experiments on expanding Bose gases and to stud- 
ies of transverse localization in nonlinear optical media. In a nonlinear medium the energy of the 
wave-packet is stored both in the kinetic and potential parts, and details of its propagation are to a 
large extent determined by the transfer from one form of energy to the other. A theory describing 
the evolution of the wave-packet has been formulated in [G. Schwiete and A. Finkel'stein, Phys. 
Rev. Lett. 104, 103904 (2010)] in terms of a nonlinear kinetic equation. In this paper, we present 
details of the derivation of the kinetic equation and of its analysis. As an important new ingredient 
we study interparticle-collisions induced by the nonlinearity and derive the corresponding collision 
integral. We restrict ourselves to the weakly nonlinear limit, for which disorder scattering is the 
dominant scattering mechanism. We find that in the special case of a white noise impurity potential 
the mean squared radius in a two-dimensional system scales linearly with t. This result has previ- 
ously been obtained in the collisionless limit, but it also holds in the presence of collisions. Finally, 
we mention different mechanisms through which the nonlinearity may influence localization of the 
expanding wave-packet. 

PACS numbers: 71.10.Ay, 71.10.Pm, 75.40.Cx 



I. INTRODUCTION 

Currently, much attention is devoted to experiments 
studying the dynamics of a wave-packet evolving in the 
presence of both random scatterers and nonlinearity. 
These experiments are inspired by the idea that one can 
visualize the phenomenon of Anderson localization. The 
propagation of a wave-packet in the presence of multiple 
scattering on a random potential has been studied using 
photonic crystals^! and also ultra-cold Bose gases con- 
fined initially inside a trarJ^HUl. The nonlinearity in the 
case of photonics is induced by the Kerr effect (the change 
in the refractive index in response to an electric field) , or 
may result from the particle-particle interactions in the 
case of cold atoms. In the optics experiments, a laser 
beam is sent into a nonlinear optical crystal with a re- 
fractive index varying randomly in the plane transversal 
to the direction of the pulse propagation. The resulting 
beam profile can be monitored on the opposite side of 
the crystal. In a second class of experiments, atoms con- 
densed initially inside a trap arc released and, during the 
subsequent expansion, are subjected to a disorder poten- 
tial. Unlike in the case of photonic crystals, in the latter 
experiments it is possible to extract information about 
the full time-evolution of the expanding wave-packets. 

Motivated by these experiments, we recently presented 
an effective theory of the propagation of a wave-packet 
(averaged over many disorder-realizations) injected in a 
disordered and nonlinear medium in two dimensions^. 
In the regimes preceding Anderson localization, or when 



it is absent, we found that the propagation of the wave- 
packet in a nonlinear disordered medium exhibits inter- 
esting dynamics related to the fact that in the presence 
of nonlinearities the energy of the wave-packet is stored 
both in the kinetic and potential parts. Then the propa- 
gation of the wave-packet is to a large extent determined 
by the transfer from one form of the energy to the other. 

The derivation of the kinetic equation presented in 
Ref. H31 was based on a classical field theory, supple- 
mented with the use of the quasiclassical approxima- 
tion, a well-known tool in the theory of nonequilibrium 
superconductivitjG^^. The corresponding functional 
can also be used as a basis for a diagrammatic pertur- 
bation theory. The relation between the different terms 
appearing in the kinetic equation and the diagrammatic 
perturbation theory was explained in Ref. 1171 Recently, 
the kinetic equation was re-derived in Ref. 18 using a di- 
agrammatic approach. In this article, we present details 
of the microscopic approach used for the derivation of 
the kinetic equation presented in Ref. Q21 We also in- 
clude an important new ingredient into the formalism, 
inter-particle collisions. As a consequence, the resulting 
kinetic equation contains an additional term, the collision 
integral. We finally discuss the relevance of the collision 
processes. 

We will assume that the time evolution of the injected 
wave-packet is governed by the nonlinear Schrodinger 
equation (NLSE), which is referred to as the Gross- 
Pitaevskii Equation (GPE) in the context of atomic Bose- 
Einstein condensates. The NLSE/GPE differs from the 
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conventional Schrodinger equation by an additional cubic 
term (we set H = 1 for the GPE): 
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= - ^ V 2 *( r , t) + u(r)*(r, t) + A|*(r, t)| 2 *( r , t). 

For negative (positive) A the nonlinear term is of the 
self-focusing (de-focusing) type. This corresponds to an 
attractive (repulsive) potential A | \E f (r, t) | 2 . The disorder 
potential u(r) is the source of randomness in the above 
equation. Starting from the NLSE/GPE, we derive a ki- 
netic equation that describes the diffusive evolution of an 
injected wave-packet in a disordered nonlinear medium. 
Since the disorder we study is static, the kinetic equa- 
tion preserves not only the integrated intensity/number 
of particles, but also the energy carried by the diffus- 
ing wave-packet. For a repulsive nonlinear term in the 
NLSE/GPE (that is typical for cold atoms), the potential 
energy stored in the medium is positive. Then, during 
the course of expansion, the potential part of the en- 
ergy is gradually converted into the kinetic part, thereby 
increasing it. For an attractive nonlinearity (typical for 
optics) , the potential energy stored in the medium is neg- 
ative, and the dynamics is richer and may, in principle, 
include a collapse^HSI] 

The NLSE used in optics is derived from the Maxwell 
equations using the so-called paraxial approximation, 22 
and thus describes the evolution of the smooth envelope 
of the electric field. The propagation direction of the 
laser beam, say the z-direction, plays the role of time 
in the NLSE. In this sense, the disorder potential which 
results from random variations of the refractive index 
is static when it is z-independent (only such a system 
is considered here). For example, the two-dimensional 
(2D) transverse evolution of a pulse is studied in a 3D 
sample 23 . In optics, the mass m has to be replaced by the 
wave vector k = uj/c, where u> is the frequency of the car- 
rier wave and c the velocity of light in the medium. The 
intensity of the beam is proportional to |\&(r, z)\ 2 . We 
will be interested in the description of the wave-packet 
when its size L = L(z) exceeds much the typical mean 
free path, Z t yp ; which in turn is much larger than the typ- 
ical wave-length Xt yp of the components constituting the 
wave packet: 



L > Z tyP > A 
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All three scales are related, of course, only to the propa- 
gation in t he dir ections transverse to z. 

The GPEP^Sis commonly used for the description of a 
large ensemble of Bose-atoms confined inside a trap. We 
are, in turn, interested in the evolution of a cloud in which 
atoms are scattered by a random potential. The usage 
of the GPE in this context is worth commenting: The 
Schrodinger equation for the field operators describing a 
many-body system, ip(r,t), can be written as 



where, under the assumption that the scattering length 
a s is the shortest length in the problem, the poten- 
tial of the particle interaction can be taken in the form 
U(r) — X5(r) (recall that for atoms A = 4TTh 2 a s /m, 
where a s is the scattering length). We will assume that 
occupation numbers n p for the relevant momenta are 
large to ensure high occupancy. In this case the opera- 
tors ip in this equation may be substituted by a complex 
valued classical field ^ (for a formal discussion of this 
point see, e.g., Ref. I26p . It is worth mentioning that in 
the case of quantum electrodynamics a similar step leads 
to the classical Maxwell equations for large photon occu- 
pation numbers. It will be important for us that the field 
\l/(r, t) should not necessarily be interpreted as a conden- 
sate wave function in order to be described by the GPE. 
The density of the cloud can be expressed as |W(r,t)| 2 . 

In addition to the condition of Eq. (1.2), throughout 
this paper it will be assumed that 
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V 2 ?A + u(r)ip + X^ijjtj), (1.3) 



where a is the inter-particle distance of atoms in the 
cloud. The former inequality corresponds to a high oc- 
cupancy of atoms which justifies the use of the classi- 
cal GPE for the description of the Bose gas. The lat- 
ter inequality means (by definition) that the gas is di- 
lute. Since we study the effects of the nonlinearity, we 
are nevertheless interested in a situation for which the 
gas is sufficiently dense in the sense that the energy per 
atom induced by the nonlinearity, which is of the order 
of A| x f'(r, t)\ 2 , is not negligible compared to the typical 
kinetic energy of the atoms constituting the cloud. 

In line with most of recent experiments on cold 
atoms/photonic crystals, we will study the den- 
sity/intensity averaged over many realizations of disor- 
der. Correspondingly, we are interested in the evolution 
of the wave-packet on length scales exceeding the typi- 
cal mean free path l typ . To obtain an averaged descrip- 
tion for the propagation of the cloud, one needs to in- 
troduce the smooth disorder averaged density, n(r,t) — 
(|*(r, t)\) dis . As a result, the nonlinearity generates a 
term of the form 2Xn(r,t)^>(r,t), i.e., it gives rise to a 
self-consistent potential $(r, t) = 2Xn(r,t). We would 
like to stress that while the density n(r, t) is smooth on 
the scale of the mean free path, the wave function ^(r, t) 
is not. Indeed, in the case we study the wave function 
varies rapidly on this scale, since the wavelength is as- 
sumed to be much smaller than the mean free path. A 
similar-looking term, 2An(r, i)^(r, i), arises in the de- 
scription of a coupled system of condensate and non- 
condensate particles, where n stands for the density of 
non-condensate pa rticle s, while \& is the smooth conden- 
sate wave function! 27 ^ 28 ! In contrast, in our description n 
is the density of the whole gas. 

The self-consistent potential is not the only effect origi- 
nating from the nonlinearity that contributes to the effec- 
tive kinetic theory of wave-packet propagation. Indeed, 
in the next order in the nonlinearity A, the so-called col- 
lision integral arises, which describes inter-particle colli- 
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sions. We will discuss this issue for atoms for which the 
meaning of collisions is more obvious. To get an idea 
about the collision rate, let us first consider the rate of 
two-body collisions in the gas of small density, for which 
the occupation numbers are small, n p <C 1. In the three- 
dimensional case, the collision rate is the inverse of the 
Maxwell-Boltzmann collision time: 1/tmb — V2n(r)crv e , 
where the atomic cross section a = Sira^. and v £ is 
the velocity of a particle with the energy e. Then, 
l/( r Af-B £ ) ~ (a s /a) 2 (A e /a), which in a dilute gas with 
small occupation numbers is a product of two small fac- 
tors (A e is the wave-length of a particle with the energy 
e). The situation changes radically for a gas with large 
occupation numbers, n p ^> 1. The smallness induced by 
the scattering length a s in the dilute gas, can be compen- 
sated by large factors n p . (The balancing between the 
smallness of the interaction amplitude and large occupa- 
tion numbers is specific for Bose-gases as compared to 
fermionic systems.) As a result, one gets for the collision 
rate 1/t co u ~ A|n 2 /e, where e is a typical kinetic energy 
of the Bose-atoms. Let us finally emphasize that while 
we used here the language appropriate for atomic gases, 
the collision rate 1 /t co u has its origin in the nonlinearity 
and as such this estimate is relevant for any system de- 
scribed by the NLSE/GPE irrespective of its microscopic 
origin. 

The kinetic equation presented in this paper is derived 
for the case when disorder is responsible for the dom- 
inant scattering mechanism, 1/r 3> 1/t co /;. To be in 
correspondence with this inequality, we will limit ourself 
to the case when the effect of nonlinearity is sufficiently 
weak so that An(r) <C e(r). 

It is worth commenting on an important byproduct of 
the interaction smallness. Under the condition An(r) <C 
e(r) we need not consider the transition to the Bogoli- 
ubov spectrum. This is because under this condition only 
a tiny fraction of the states with the smallest energies is 
influenced by the off-diagonal components in the Bogoli- 
ubov Hamiltonian. For the majority of the particles the 
off-diagonal components of the Bogoliubov Hamiltonian 
can safely be ignored. 

As it was already mentioned, when treating disor- 
der we assume that the mean free path is much larger 
than the typical wavelength \ typ of the components con- 
stituting the wave-packet. Throughout this paper, we 
use the model of a delta-correlated Gaussian disorder 
potential, characterized by (u(r)u(r')) = 7<S(r — r'). 
This model is appropriate if scattering occurs on quan- 
tum impurities, for which the range of the potential is 
much smaller than the wavelength Xt yp - For the delta- 
correlated disorder potential, the density of states deter- 
mines the frequency-dependence of the scattering rate, 
1/t{e) = 2-Kv{e)^i cx e( rf ~ 2 )/ 2 and of the diffusion coef- 
ficient D(e) — 2er(e)/mrf cx e 2 ~ d / 2 . In particular, the 
scattering rate in d — 2 is energy-independent. Both 
in optics experiments and in experiments on Bose gases, 
one often uses speckles to realize the disorder potential. 
The speckle potential has a finite correlation length. If 



the wave-length \ typ is much larger than the correlation 
length, the model of the delta-correlated disorder poten- 
tial remains a good approximation. If the wavelength is 
sufficiently short to resolve the finite correlation length, 
however, one needs to be more cautious. Unlike for the 
short range scatterers, the typical time for the random- 
ization of the momentum direction, i.e., the transport 
scattering time, no longer coincides with the single parti- 
cle scattering time, which is determined by the imag- 
inary part of the self-energy in the disorder averaged 
Green's function. The transport scattering time r tr (e) 
acquires a frequency dependence that differs from the 
one for short range scatters stated above. The same is 
true for the diffusion coefficient, since it depends on r tr 
as D = 2er tr (e) /md. The expression for r tr appropriate 
for a speckle potential can be found in the literature, e.g., 
in Refs. 1291 As concerns the nonlinear diffusion equation 
derived in this manuscript, it can be expected that the 
only change that needs to be introduced when dealing 
with a speckle potential is the replacement t — > r tr in 
the final form of the equations, which already contains a 
energy-dependent diffusion coefficient. 

The paper is organized as follows. In Sec. [TT] we pro- 
ceed directly to the discussion of the nonlinear kinetic 
equation. Those readers, who are not interested in the 
technical details of the derivation of the kinetic equation 
based on the quasiclassical approximation, find the most 
important information in Sec. [n] as well as in the Con- 
clusion. First, we discuss the equation in the collisionless 
regime in Sec. |II A| Although most of the material of 
Sec. |II A| has already been presented in Ref. [I3J we in- 
clude it here in order to make the paper self-contained. 
In the second part, Sec. |IIB| we add the effect of colli- 
sions. It turns out that the interparticle collisions impose 
certain constraints on the range of validity of the derived 
equations. The main result of this paper is formulated 
here: In two spatial dimensions, the mean squared radius 
of the wave-packet grows linearly in time. This result is 
not affected by inter-particle collisions. 

In Sec. |III| we introduce the field theory approach that 
is the main tool for our investigations. The basic idea is 
to write a functional integral expression for the time evo- 
lution of the observable in question. (Our aim here is to 
describe the evolution of the density/intensity n = \^\ 2 . 
The wave function at the initial time 'J'o is assumed to 
be known.) Typically, this kind of approach is used when 
studying Langevin-type equations including a noise term 
with a given correlation function. In the problem un- 
der study in this paper, no noise is considered. Instead, 
we use an analogous construction, and then average over 
disorder configurations. The resulting theory closely re- 
sembles the structure one encounters in Keldysh field the- 
ories, where Green's functions can be transformed to a 
block-triangular form. Retarded and advanced Green's 
functions are supplemented by a third type of Green's 
function that contains information about the distribu- 
tion function n(r, t, e), which we are interested in. 



In Sec. IV the averaging over the disorder potential is 
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performed, i.e., we provide a description of the evolution 
of a wave-packet averaged over many disorder configura- 
tions (realizations). First, the theory of the wave-packet 
in the absence of the nonlinearity is discussed. Here, 
we make contact with Ref. [301 and 1311 where the expan- 
sion of a Bose-condensate over a disorder potential was 
studied starting from a later stage of the time-evolution 
when the nonlinearity may already be neglected. After- 
wards, the nonlinear problem is considered. We start 
this discussion with a diagrammatic analysis (in two di- 
mensions) before deriving the kinetic equation using the 
method of quasiclassical Green's functions. Here, we pro- 
ceed in close analogy with the theory of nonhomoge neo us 
superconductivity^^. The main result of Sec. IV 



Note that the diffusion term contains a sort of the co- 
variant derivative: 



is 



given by Eq. (2.5 1, which is a classical nonlinear diffu- 



sion equation in the collisionless regime. The equation 
was first presented and analyzed in Ref. [13] for a two- 
dimensional system. Discussion of two dimensions was 
of special interest for us, because for weak disorder there 
is an exponentially large diffusive regime before the An- 
derson localization takes place. After our work^, the 
equation Eq. (2.51 was re-derived and generalized for ar- 



bitrary dimensions in Ref. 1181 using the diagrammatic 
technique. It was noted that for a generalization to di- 
mensions d 2 a new term in the kinetic equation is 
required in order to account for the non-constancy of the 
density of states. In Sec. |IV| the equation is obtained for 
arbitrary dimensions d = 2, 3 including the additional 
term found in Ref. HU 

In Sec. [V] we derive the collision integral in the kinetic 
equation originating from the NLSE/GPE. We provide a 
diagrammatic interpretation of the different terms con- 
tributing to the collision integral. Finally, we conclude in 
Sec. |VI| with a discussion of the results. In particular, we 
comment on the role of the nonlinearity in the context of 
localization. 



II. DISCUSSION OF THE NONLINEAR 
KINETIC EQUATION 

A. The kinetic equation in the collisionless regime 

We start from the nonlinear kinetic equation deter- 
mining the density evolution in the diffusive regime. The 
argument e in this equation has the physical meaning of 
the kinetic energy, e(r, t) — e — #(r, t), while $(r, i) is a 
self-consistent potential. Correspondingly, the diffusion 
coefficient is Dg = v^Tg/d. Then, the equation for the 
distribution function looks as follows: 

d t n(r, t, E ) - V(£>eV r n(r, t, e)) 

+<9 f t?(r, t)d e n(r, t, e) = 6(t)2nv(e)F{i, r). (2.1) 

This equation should be supplemented with the self- 
consistency relation for the potential i?(r, t) = 2An(r, t), 
where 

n(r,t) = I J n(r,t,s). (2.2) 



V r = V- W(r,i)r e ~, 



(2.3) 



where Fg = — d e lnz/(e). The term on the right-hand side 



of Eq. ( 2.1 1 specifies the injection of the wave-packet and 



initial evolution up to times of the order of the scattering 
time t. Namely, 



F(e,r) = J 



|^F(p,q)e«" r 2ri(e^ p ), (2.4) 



and F(p, q) = ^ (p + q/2)*g(p - q/2) is determined by 
the initial wave function ^q. Further, e p = p 2 /(2m) is 
the kinetic energy. 

Despite its apparent simplicity, Eq. ( |2.1[ ) is a rather 
complicated nonlinear integro-differential equation. The 
diagrammatic interpretation of the different terms ap- 
pearing in this equation is provided in Sec. |IVB| for the 
two-dimensional case. The main new ingredient for d =/= 2 
is the non-constant density of states, v (i). (Note that T 
vanishes in two spatial dimensions when the density of 
states is constant. In three dimensions, however, T = 
is finite.) Since the density of states enters with the ar- 
gument e = e — $(r, t), the scattering rate acquires an 
explicit dependence on This eventually leads to a mod- 
ification of the diffusion term in Eq. ( |2.1[ ) by substituting 
V — > Vr , which was first noticed in Ref. (T5J 

The underlying physics of the nonlinear diffusion equa- 
tion Eq. (2.1) was discussed in Ref. [13l The equa- 



tion describes diffusion of a particle with total en- 
ergy e on the background of a smoothly varying po- 
tential #. Correspondingly, the kinetic energy e p = 
£ — $ varies locally in space and time. One may no- 
tice that in the NLSE/GPE a purely time dependent 
potential may be removed by a gauge transformation 

#(r,t) -> $(r,t)exp(-iJ* o dt'V(t'f), that leaves the 

density | ^ (r, £) | 2 unchanged. On the level of the dis- 
cussed equation, this point becomes obvious when writ- 
ing the distribution function as a function of the ki- 
netic energy instead of the total one, n(r, e, t) = h(r, e + 
#(r, t),t). Expressed in the new coordinates the equation 
reads 



d t n(r,e,t) 

- V - W(r, i)d £ 

= 8{t) F(e,r), 



D r 



V0(r,t)0 e Jn(r,e,t) 

(2.5) 



where now Vr = V — V$(r, t)T e . One can see explicitly 
that a purely time-dependent potential drops from the 
equation since i?(r, £) enters only in combination with a 
spatial derivative, as Vt?(r, t). In Eq. (2.5), the diffu- 
sion coefficient D(e) = 2£t(e) /md depends explicitly on 
e, but also implicitly through t(e). Within our model of 
a delta-correlated impurity potential, the elastic scatter- 
ing rate l/r(e) acquires a frequency dependence through 
v(e). The form of the equation suggests, however, that 
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it will also hold in the case of impurity potentials with 
a finite correlation length, when t(s) should be replaced 
by the transport scattering time T tr (e). 

It seems clear that a closed form solution of the non- 
linear equation for arbitrary initial conditions cannot be 
found. In order to make progress we will rely on the use of 
conservation laws. The GPE describes a system in which 
the total particle number (or intensity in the case of the 
NLSE) and the total energy are conserved. The total 
momentum is not conserved, since the disorder potential 
breaks translational invariance. It is important to check 
that our approximations are consistent with the conserva- 
tion laws, namely that energy and number conservation 
are still encoded in the nonlinear diffusion equation (2.5 1. 



Let us start with the number conservation. For that, 



we integrate Eq. ( 2.5 1 in e and obtain the continuity equa- 
tion in the form dtn\r, t) +Vj(r, t) — 5(t)n(r,t). The role 
of the right hand side is merely to determine the bound- 
ary condition at the initial time t = 0. The expression 
for the current is 

r de 

jM) = J ^j( r :M) (2.6) 

j(r,e,i) = -£> e [Vr-W0 e ]n(r,£,t) (2.7) 

Next we turn to energy conservation. Here, the con- 
tinuity equation, 9 t ps(r, t) — — Vj^, takes the following 
form: 



p E (r, t) = e(r,t) + An 2 (r,i) 
cfe 
2^ 



3E(r,t) 



(2.8) 

(e + 0)j(r,t >e ) (2.9) 



where e = / (de/2n) en(r, t, e) can be interpreted as the 
average kinetic energy. In particular, we may conclude 
that the total energy 



Etnt = 



J dr (e + An 2 ) 



(2-10) 



is conserved. The total energy is conserved for our prob- 
lem, because impurity scattering is elastic and we con- 
sider a closed system. The conservation of energy is a 
known property of NLSE/GPE from which we started. 
The derivation based on the kinetic equation, which we 
presented here, can be regarded as a check of the validity 
of our approach. 

Remarkably, as we have observed in Ref. [T31 for two 
spatial dimensions when T = 0, and if the scattering time 
is frequency-independent, the conservation laws com- 
pletely determine the time evolution of the mean radius 
squared of the wave-packet, (r 2 ) = J dr r 2 n(r,t)/N. 
Indeed, in 2d the expression for the current j(r, t) can be 
simplified and the continuity equation takes the form 

d t n(r, t)-— V 2 (e + An 2 (r, t)) = 5{t)n{v, t). (2.11) 
m 



Now multiplying Eq. (2.11) by r and subsequently inte- 



grating in r one obtains that 

d t (r 2 ) = AD e 



where Etot = E to t/N . The linear dependence of the mean 
square radius on time during the evolution is guarded by 
energy conservation. This is a rather non-trivial result; 
the rate of expansion is proportional not to D^, as one 
may naively expect, but to DE tot . The reason is that 
the rate of expansion combines the effect of diffusion and 
propagation in the field of the force induced by the self- 
consistent potential. This is one of the central results 
of our previous papei^, unfortunately, in higher dimen- 
sions it seizes to be valid due to the non-constancy of the 
density of states. 

It remains to discuss general features of wave-packet 
dynamics in the repulsive and the attractive case. When 
the potential energy related to the nonlinearity is con- 
verted into kinetic energy, the total kinetic energy in- 
creases in the repulsive case and decreases in the attrac- 
tive case. Correspondingly, during the course of the ex- 
pansion localization effects can be expected to be weak- 
ened for repulsive nonlinearity and enhanced for attrac- 
tive nonlinearity. In particular, for an attractive non- 
linearity the slowing down and eventual localization of 
the injected pulse (not considered here) occurs at smaller 
distances than in the linear case as observed in the 
experiment!^. As it was indicated in Ref. [T31 if a part of 
the cloud lags behind, this fragment may have a strong 
tendency to localize. One may expect that this kind of 
localized fragment generically remains from an expand- 
ing cloud when the nonlinearity is attractive. To check 
this point, it would be desirable to analyze data with 
respect to the intensity/number of particles of the re- 
maining localized part of the cloud and, if possible, the 
energy concentrated in this part as compared to that in 
the initial cloud. 



B. The role of collisions 

The nonlinear term in the NLSE/GPE gives rise to a 
collision integral in the diffusive kinetic equation, which 
is proportional to A 2 . The full kinetic equation including 
intcrparticle collisions takes the form 



d t n(r,t,e) 

- V - Vtf(r,t)<9 £ 



V r -Vt?(r,i)0 E 



= 5{t) F(r,e) + 2 7 rzy( £ )/ coll (r,i,e) 



n(r,t,e) 
(2.13) 



with 



rColl 



(r,f,e) 



— in X 2 (2n) d J dndn2dn3<iri4 J dezdezde^ 

xv(£ 2 )v{e 3 )v(e i )5{e + e 2 - e 3 - e 4 ) 
xd(p e + p E2 -p £3 -Pej([n' e +n' E2 \n' e3 n' ei 



-n £ n S2 [n es +n £4 . 



(2.14) 



(2-12) 



where 2irv(e)n' E (r, t) = n(r,t,e), p E = p E n, p £i = p ei n iy 
and the integration goes over positive frequencies only. 
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To conclude, we get a standard collision term of two par- 
ticles in the limit of large occupation numbers n' E . 3> 1. 
The left-hand side of the kinetic equation takes into con- 
sideration that the distribution function of states partic- 
ipating in the collision are determined by the diffusive 
propagation in the disordered and nonlinear medium. 

The collision integral contains two terms describing 
the "in"- and "out" -collision channels. To estimate the 
scattering rate l/r co u, let us focus on the "out"- term, 
which is given by the last term in the expression for I co ih 
Eq. (2.14), and is proportional to n' £ . We will write it 
as n e /? 'coll- Recall that the typical kinetic energy per 
particle at point r is denoted as e(r). For a conservative 
estimate of the scattering rate, let us consider an energy 
e ~ e(r); in this case the kinematic constraints induced 
by the momentum and energy conservation in the col- 
lision integral are minimal. Since one has to integrate 
two distribution functions over energies, this ultimately 
yields a factor n 2 (r). As a result one gets 



1 

Tcoll 



A 



2 n 2 (r) 
e(r) ■ 



(2.15) 



It is clear from this estimate that in order to use the 
language of the kinetic equation with well defined dis- 
tribution function n(r, t, e), one has to be limited to 
the case when Xn(r,t) <C s(r). Under this condition, 
1/t co u < e(r). 

Still, there remains a question about a comparison 
of the rate of inter-particle collisions with elastic scat- 
tering caused by disorder, i.e., 1/t co u versus 1/t. In 
this paper we limit ourself to the case of rare collisions, 
1 / 1 T C oU "C 1/t, i.e., we assume that elastic scattering 
events occur more frequently than inter-particle colli- 
sions. This condition is more restrictive than the con- 
dition e(r) ^> 1/ T co u discussed above. 

The collisions, naturally, change the dynamics of the 
propagation. As long as the kinetic eq uation in the de- 
rived form holds, however, the result (2.12) about the 



rate of the expansion of the wave-packet remains valid 
even in spite of the inter-particle collisions. This is be- 
cause (i) the collision integral is local and as such does 
not change the mean radius squared of the wave-packet, 
( r2 ) dis - Furthermore, (ii) in two spatial dimensions, the 
rate of expansion depends only on the total energy E to t, 
which is not altered by collisions and it docs not depend 
on the energy dependence of the distribution function, 
which is controlled by the collision integral. 

Finally, we would like to note that while the rate of 
"delivery" of colliding particles was controlled by diffu- 
sion, we did not consider the modifications of the col- 
lision integral by disorder. It is is very different from 
what happens in disordered conductors at low temper- 
atures, T <C 1/r <C £f- The reason is that the ki- 
netics of the classical particles, not constrained by the 
existence of the Fermi-surface, is similar to the case for 
which 1/t<T~£f with ef ~ £, where Ef is the Fermi 
energy. Then, modification of the collision integral by 
disorder leads to a smallness l/re(r) without gaining a 



large factor l/(rT), as it was in the case of conductors 
at low temperature. 



III. BASIC FORMALISM 

In this section we introduce the field theory approach 
that is the tool for our investigations. Our aim is to 
describe the evolution of the density (intensity) n = \^\ 2 , 
averaged over disorder configurations. The wave function 
at the initial time \l/o is assumed to be known. 

Formally, the problem bears a certain similarity with 
the description of critical dynamics near a phase tran- 
sition, or, more generally, the study of Langevin-type 
equations with the help of field theory approaches. The 
formalism we are alluding to here i s often called Martin- 
Siggia-Rose (MSR) formalis m 1 28 1 321 ^ an d finds applica- 
tions in many different branches of physics. The basic 
idea is to write a functional integral expression for the 
time evolution of the observable in question. With the 
help of a delta-function entering the integral, the wave 
function is fixed to coincide with the solution of the un- 
derlying equation. By introducing an additional field 
variable and thereby doubling the degrees of freedom, 
one may write the delta-function with the help of an in- 
tegral over an exponentiated action. 

Typically, this kind of approach is used when studying 
dynamical problems, for which the original equation con- 
tains a noise term with known correlation function. One 
may then average over the noise, and study the resulting 
functional with field theoretical methods like perturba- 
tion theory, the renormalization group, or by analyzing 
instantonic configurations. In the problem under study 
in this paper, no noise is present. Instead, we use an anal- 
ogous construction, and then average over disorder con- 
figurations. With a proper regularization, vacuum loops 
are absent right from the beginning and this is why the 
dynamical approach is particularly useful for the prob- 
lem of quenched disorder, as was already noted long time 
bacfeP. 

The resulting field theory indeed closely resembles 
the structure one encounters in Keldysh field theories, 
where Green's functions can be transformed to a block- 
triangular form. Retarded and advanced Green's func- 
tions are supplemented by a third type of Green's func- 
tion that contains information about level population. 

For a Bose-Einstein condensate, one can obtain the 
Gross-Pitaevskii equation as a mean field equation for 
the full quantum many-body problem. As one might ex- 
pect from this observation, a connection exists between 
Keldysh-type field theories for the quantum problem, and 
the MSR-type approach. Indeed, in the Keldysh ap- 
proach, two distinct types of interaction vertices exist, 
they are sometimes referred to as quantum and classical 
vertices^. By disregarding the quantum vertices, while 
retaining the classical ones, one recovers a representation 
of the functional delta-function, that fixes the evolution 
of the (classical) fields to obey the classical equation of 
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motion, in this case the Gross-Pitaevskii equation. This 
approach additionally allows to consider correlations in 
the initial density matrix, and one can obtain, for exam- 
ple, the so-called truncated Wigner approximation, as ex- 
plained in more detail in Ref. 1351 In optics, the nonlinear 
Schrodinger equation emerges as a result of the paraxial 
approximation applied to the Helmholtz-equatiorP 2 ! and 
has thus a different microscopic origin. This is the reason 
why we do not explicitly use the (microscopic) Keldysh 
approach as a starting point in this paper. 



A. Action 

Our starting point is the Gross-Pitaevskii equation in 
the form given in Eq. (1.2 1. This equation describes the 



time evolution of the macroscopic wave- function ty(r,t) 
in the presence of an external potential it(r). The total 
density |\P(r,t)| 2 is conserved in time and we use the 
normalization J dr] 1 3 r (r, t)| 2 = N, where N is the total 
number of atoms in the gas. The quantity of our interest 
is the disorder averaged density 



n(r,i) = (|*(r,t)| 2 ) 



(3.1) 



Disorder averaging (. . . ) dis is performed with the help of 
the Gaussian probability distribution 



V(u) = Afexp 



■- / dr u(r)W _1 ( 



r>(r') ,(3.2) 



where J\f provides the normalization, so that 
/ Du V{u) = 1. This definition implies that (u(r)) dis = 
and (u(r)u(r')) = W(r — r'). In this paper, we consider 
the specific case of a delta-correlated (white noise) 
potential, for which W(r — r') = j8(r — r'). In two 
dimensions the density of states v is constant, and one 
can identify 7 = 1/2i:vt, where t is the scattering time. 

We first note that the unaveraged density can be rep- 
resented as the following functional average: 



n(r,t) = / D(^,r)D(v,V*) l*M)| 



2 AS 



where we introduced the complex fields 77, rj* and ip, ip* 
The action S is given by 



5' 



drdr'dtdt' 



r?*(r,i) 



V>(r,t) 
V(r,t) 



+i / dr [Tj(r,0)*5(r)-tj*(r,0)* (r)], 



(3.4) 



where the inverse matrix Green's function g 1 has the 
structure 



The retarded and advanced Green's functions Qr/a fulfill 
the equation 

(idt + ^ - «(r) - A|V(r, t)\ 2 ^j g R/A (r, r', t, t') 

= 5(r-r')5(t-t') (3.6) 

with standard boundary conditions. Indeed, upon inte- 
gration in the auxiliary fields rj(r,t), rj*(r,t) one obtains 
a functional delta function that fixes the fields \P(r, t) 
and fy*(r,t) to obey the Gross-Pitaevskii equation and 
its complex conjugate, respectively. The last part of the 
action involving the fields ^0 and ^0 fixes the bound- 
ary conditions at the initial time, ^(r, to) = v I / o(r) and 
v I'*(r, to) = x I , o( r )- We see that the formalism involves a 
doubling of the degrees of freedom, similar to the Keldysh 
or closed-time-path approaches for quantum systems^, 
where two fields are introduced on forward and backward 
time-contours. We repeat that with a proper regulariza- 
tion vacuum loops are absent. For a more detailed ac- 
count of the construction of the classical funct ional an d 
the appropriate regularization we refer to Refs! 28 ! 36 * 37 ! 

In order to lighten the notation, we find it conve- 
nient to introduce the field doublets <f» = (tp, rf) T and 
(j) = (fta x = (r]*,tp*), so that 

n(r,t) - (tr(<7_ [0(r,t) ® 0(r,t)])) (3.7) 

Pauli matrices <jj act in the space of the the fields ip 
and 77, and <r_ = {a x — ia y )/2. The averaging (...) — 
J D((f>, (fy) (...) exp(iS') is performed with respect to the 
action S, which we write in terms of (f> and <f> and split 
into several parts, 



S — So + S s 



u dis 



(3.8) 



The term So alone describes the free propagation of fields 
4> in the absence of interactions and impurities. The 
source term S s contains information about the initial 
conditions, for convenience we choose to = from now 
on. The disorder potential and the nonlinear term in 
the Gross-Pitaevskii equation give rise to S' dis and Si nt , 



(3.3) respectively. 



So = / drdv'dtdt' ct>(r,t)g- 1 (r-r',t-t')cj)(r',^ 1 9) 



i I dr (0 o (r)^(r, 0) - 0(r, 0)^o(r)) , (3.10) 
drdt^(r,t)u(r)cj){r,t), (3.11) 
-A / drdt 0(r,t)0(r,t) ~cj>(r,t)<T_<f>(r,t). (3.12) 



S s 
qi 

u dis 
Sint 

Here, the 2x2 matrix Green's function 



g-i _/ 5a 1 



(3.5) 



a - (^ ° 



(3.13) 
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is composed of the retarded and advanced Green's func- 
tions gy A (p,s) = (e — e p ± i<5) _1 , where e p = p 2 /2m. 
For the initial condition, we introduced 

o (r) = (* o (r),O) T , o (r) = (O,*5(r)). (3.14) 

It is an important consequence of the structure of the 
theory, that the Green's function G = — i (0 0) h as a 
triangular structure, where in accord with Eq. ( |3.13 ) the 
11 and 22 elements are retarded and advanced Green's 
functions. These Green's functions contain information 
about the spectrum, while the off-diagonal (12) element 
contains information about the occupation, in analogy 
to the Keldysh approach. Importantly, the 21-element is 
equal to zero. 



B. Diagrammatic representation 

We start with an elementary discussion of the structure 
of the perturbation theory. We will draw diagrams in 
such a way that time runs from left to right. Retarded 
and advanced Green's functions are depicted in Fig. [T] 



'/ 

•- 



G R 
-> — 



1> 



V 
•- 



G A 



FIG. 1: The retarded (G R ) and advanced (G A ) Green's func- 
tions. The time arrow runs from left to right. 

The close similarity to a Keldsyh field theory has al- 
ready been stressed above. The main difference com- 
pared to a full quantum theory of interacting bosons in 
the Keldysh approach is that out of the two types of ver- 
tices depicted in Fig. [2] only one is realized. Namely, only 
the so-called classical vertices, shown on the left hand 
side of Fig. [2] appear in the theory considered here, while 
the so-called quantum vertices, shown on the right hand 
side, are absent (see the related discussion in Ref. I3"5"j) . 
This has important consequences. It immediately implies 
that the interaction vertices related to the nonlinearity 
have the structure shown in Fig. [3j This structure, in 
turn, implies that there are no closed loops in this rep- 
resentation. In order to draw more complex diagrams 
in a convenient way, we will often depict the interaction 
vertices with an additional wiggly line (as, for example, 
in Fig. [1] below), but one should keep in mind that the 
interaction is in fact local in space and instantaneous. 

In order to further elucidate the structure of the per- 
turbation theory, we study the expression for the density 
evolution. The disorder averaging is postponed until the 
next section, in this section all Green's functions are un- 
averaged and explicitly depend on the disorder potential. 
First, we introduce two real Hubbard-Stratonovich fields 
fid and "d q , which we assemble into the following matrix: 



d 



$ q #cl 



With the help of this matrix, the interaction can be rep- 
resented as 



exp (iS int ) = (exp(iS#))0 , 
where we introduced the notation 



(3.16) 

dvdt 0(r,i)tf(r,i)0(r,i), (3.17) 



and (. . . )ft symbolizes the the following averaging proce- 
dure 

{...)+ = J D-d (...)eA/ drdtflT ( r . t )^''(^).(3.i8) 

In this equation, •& = (i? 9 , $ c ;) and J\f is a normalization 
constant which we will suppress from now on. 



Formula (3.18) implies that fields i} q and f) c i couple 
to each other, but not among themselves. The field "d c i 
enters S$ like a classical potential. The quantum compo- 
nent dq couples retarded and advanced Green's functions 
in a specific way. Taken together, these observation imply 
that all possible diagrams have the structure indicated in 
Fig. [4] It is also instructive to further integrate in </>. The 
result is 

n(r,t) 

J drxdr 2 *S(ri)G^ i (r 1 ,r;0,t)Gf ci (r,r 2 ;i,0)* (r 2 ) 

Xe ifdr 3 dr 4 Vo(r 3 )[G£ cl »d q »G$J(r 3 ,a,r 4 ,a)y (r A )\ , g ^ 

We used the triangular structure of G in order to obtain 
this result. The filled circle • symbolizes a convolution 
in space and time. The retarded and advanced Green's 
function in the presence of the classical field G$ cl fulfill 
the differential equation 

(idt-H-fidfatj) G^f(r,r',t,l/) 

= S(r-r')5(t-t'). (3.20) 

and H = -V 2 /2m + u(r). Before averaging in §, the pre- 
exponential factor describes the evolution of the density 
on the background of an external classical potential . 
The exponential contains a similar structure: Each term 
in the expansion of the exponential symbolizes the evo- 
lution of the density up to a certain point. From a for- 
mal perspective, integration in "8 q introduces a functional 
delta function, that fixes i? c / to equal the density. 




77* V* 





n ti 




(3.15) 



FIG. 2: In a Keldysh many-body approach to interacting 
bosons two classes of vertices appear, the classical vertices 
shown on the left and the quantum vertices shown on the 
right. In the MSR-type approach used in this paper only the 
classical vertices are present. 
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Averaging is performed with respect to the action S after 
the disorder averaging, i.e., self-cosistently. This implies 
that, generally speaking, the disorder part of the self- 
energy also implicitly depends on the interaction (namely 
via the Green's function —% (cjxj))). 



FIG. 3: Upon averaging with respect to the fields ip and 77, 
the (classical) interaction vertices in our approach give rise to 
the two sub-diagrams shown above. 




FIG. 4: General structure of the perturbation theory: The 
density evolution is represented by the infinite sum of all di- 
agrams of the type displayed in this figure. 



IV. TIME EVOLUTION IN A DISORDERED 
MEDIUM 

Initially, the disorder potential u(r) enters the action in 



the form S' 



(lis 



J drdt <fi(r,t)u(r)(f)(r,t). Disorder av- 



eraging with respect to the probability distribution (3.2) 
introduces an effective interaction of the fields 



Sdis 



7 / dvdt 1 dt 2 0(r,*i)</>(r,£i)0(r,t 2 )</>(r,t2 



(4.1) 



This effective interaction is local in space, but non-local 
in time. 

It is usually not possible to take into account disorder 
effects exactly and one needs to employ approximation 
schemes. Disorder averaging introdu ces a quartic term 
in the action S, namely Sdis of Eq. (4.1 1. Here we will 



treat this term in the self-consistent Born approxima- 
tion (SCBA), which relies on the weak disorder condition 
et ^ 1, where e is the characteristic scale for the kinetic 
energy in the problem. 



The SCBA consists in replacing Sdis given in Eq. (4.1 ) 



by 



S dis (4.2) 
17 / drdt l dt 2 <j)(r, h) (</>(r, *i)0(r, t 2 )) </>(r, t 2 ). 



The average can be taken in two equivalent ways, which 
explains the additional factor of 2 compared to Eq. (4.1 1. 



A. Noninteracting theory 

This section contains an elementary discussion of the 
theory for the density evolution in the noninteracting 
case A = 0. It serves as a preparation for the discussion 
of the interacting model. Furthermore, we use the oppor- 
tunity to introduce our notation and to stress the most 
important differences to the calculation of the density- 
density correlation function in disordered electron sys- 
tems. 

In this case one obtains 

n(r,t) = j dr l dr 2 {%(v 1 )G^(v 1 ,T;Q,t) (4.3) 
G£(r,r 2 ;t,0)* (r 2 )} dis 

where Ud t - H\ G^ M (r, r', t, t') = S(r - r')S(t - f). 

In the SCBA, the disorder averaged Green's function 
is given by 



where 



2iru(e)j 



(4.4) 



(4.5) 



is the scattering rate and v(e) is the density of states. 
This result is obtained as follows. For A = 0, the defining 
relation for the disorder part of the self-energy in the 
SCBA is 



(dp)- 



1 



~*dis 



(4.6) 



Here and in the following we use the notation (dp) 
d p/(2%) d . The scattering time t £ is defined as 



3[£^( £ )] = T-l/(2r E ) 



(4.7) 



Upon introducing the variable £ p = e p — e the integration 
measure transforms as J (dp) = J^ £ d^ v v(e + £ p ), where 
the trivial angular averaging has already been performed. 
Focusing on the imaginary part of the self-energy first, 
one may extend the lower limit of the integration in £ p to 
—00 in the weak disorder limit, e — 3?[S^ s (e)] >• l/r E . At 
the same time, this step regularizes the integral for the 
real part of the self-energy. The integrand for the imag- 
inary part of the self-energy is strongly peaked around 
£ p = and one may replace v(e + £ p ) ~ v(e) and take 
the density of states out of the integral. The remaining 
integral is easily performed and the result is 



(4.8) 
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FIG. 5: Diagrammatic representation of the diffusion process 
in the absence of the nonlinearity. 



in agreement with (|4.5j) and (4.7). 

As is well knowrP^in the leading approximation in 
1/gT , one should not only replace Go by G in formula 
(4.3) for the density, but sum the whole set of diagrams 
with non-crossing impurity lines as shown in Fig. [5] Ef- 
fectively, this amounts to summing a geometric series. 
This procedure leads to the expression 



n(q,w)= / (dp)(de) * (p+)* (P-) 

oo 

xG fl (p+, £+)G A (p_, e_) £"(q, w) (4.9) 



n=0 



where 



L e (q,w) - 7 J(dpi)G R ( Pl+ ,e + )G A ( Pl -,e-) 

We use the notation (de) = de/(2n) for frequency inte- 
grals, p± = p±q/2 and e± = e±u>/2. The expression is 
quite similar to the familiar density-density correlation 
function in electronic systems. Note, however, that in 
the latter case the frequency integration is restricted to 
a small interval around the Fermi surface of order of the 
temperature by the presence of a distribution function. 
In contrast, here the momentum integration is restricted 
by the initial wave functions an d ^q- The frequency 
integration, on the other hand, is a priori not limited. 

Let us assume that the inequalities er £ 3> 1, ujt £ <C 1, 
ql £ <C 1 are fulfilled (diffusion approximation), where 
l £ = v £ t £ is the mean free path, v £ = p £ /m and p £ — 
\Jlms are the velocity and the momentum at energy e. 
In this case we can calculate the sum approximately by 
using the expansion 



L e (q,cj) « 1 + jwt £ 



Z £ V/2. 



(4.10) 



It will be useful to introduce a frequency dependent dif- 
fusion constant as D £ — v £ r £ /d in dimension d. After 
performing the sum in the equation for the density we 
obtain 



n(q,w) = / (dp)(de) * (p+)*o(P-) 



1 



(4.11) 



xG"( P+ ,e + )G"(p_,e_)-P £ (q, W ), 

T £ 

where the energy dependent diffuson is 

V £ (q,u) = (D £ q 2 -icu)- 1 . (4.12) 



The next step is to integrate in e, where one encounters 
the following integral 



(de)- 



1 



1 



1 



-p+ 



2r e 



z> E (q,w) 



2t> 



(4.13) 



For e p+ ~ £ p 3> 1/r the most important e are of the 
order of s p and we can perform the integral with the help 
of the residue theorem considering the poles originating 
from the Green's functions only, thereby effectively re- 
placing D £ by D £ . A distinction between e p+ and e p _ 
in the argument of the diffusion coefficient would be be- 
yond the accuracy of our approach. The result is 



7i(q,w) 



J (dp)F(p,q)2? 8p (q,w), 



where we introduced the notation 

F(p, q) = * (P + q/2)*5(P - q/2). 



(4.14) 



(4.15) 



It is clear from the previous arguments that the approach 
is valid as long as £ p r e S> 1. Typical momenta p are 
controlled by the initial wave-function ty. For the aver- 
aged density as a function of coordinates and time we 
find the expression 



n{r,t) 



(4.16) 



/ (dp) 4^6 / ^ e ^ ri)2/i4D ' pt)F (P^)- 

For |ri| <C |r|, i.e. for distances |r| exceeding by far the 
extension of the initial wave-packet, we may neglect ri 
in the exponent and obtain 



*(r,t) = e(t) [(dp) 



\Mp) 



(4.17) 



This expression was presented in Ref. 30. 

In the calculation described in this section, the 
frequency-integration was performed before the momen- 
tum integration in p (see Eq. 4.13). Relevant momenta 



in the integral of Eq. (4.14) are determined by F(p, q), 
which encodes the information contained in the initial 
wave-function >f '. For a generalization to the interacting 
case, it will be more useful to perform the integration in p 
before the integration in e. In order to achieve this goal, 
we introduce the distribution function / in the following 
way 

/(r,fi,f a ) = 7 J dr 3 dr 4 G R {r 1 -r 3 ,t 1 ) 

x*(r 3 )<T (r 4 )G A (r 4 - r u -t 2 ). (4.18) 

It describes the initial section of the diffusion ladder, 
compare Fig. [5] With the help of this definit ion one can 
write 



(4.19) 



n(r,t) = / (de)n(r,e, t), 



11 



where the energy resolved density is 



n(r, E ,t) =2ttv(e) [ 2? e (r-n,t-ii)/(ri,e,t)(4.20) 



and 



f(r,e,t) = J d(At)f(r,t + At/2, t- At/2) e ieA *.(4.21) 

We can make contact with the previous results of this sec- 
tion by noting that for times t 3> t one can approximate 
(see Appendix |A| 



27ri/(e)/(r,e,t)«*(i)F(e,r), 



(4.22) 



F(e, r) = / (dp) F(p, r) (2tt)5(s - e p ). (4.23) 



where 



This concludes our discussion of the non-interacting the- 
ory. 

As for the electronic systems, it is most convenient 
to formulate the microscopic theory with the help of 
a frequency dependent distribution function, since mo- 
mentum is not conserved during the scattering process. 
At the same time, the initial distribution is determined 
by the momentum dependence of the wave function, 
Eq. (4.15). In the quasiparticle approximation, one can 
translate between the two representations, Eq. (4.231. 



The specifics of the given problem in comparison with 
diffusion in a degenerate electronic system is that the 
dependence of the diffusion coefficient needs to be kept 
explicitly. Each particle at a given energy diffuses with 
its own diffusion coefficient and the total density is ob- 
tained through a convolution with the distribution func- 
tion, Eqs. (4.16) and (4.20). The fact that the energy dis- 



tribution may be broad has the important consequence 
that the density may differ considerably from the form 
n(r,t) oc exp(— cr 2 /t) (with a constant c), which holds 
for diffusion at a fixed energy. To illustrate this impor- 
tant point, we briefly discuss an example first introduced 
in Ref . M 



1. Gaussian initial distribution 

As an instructive example, one can easily calculate the 
asymptotic distribution for the initial conditional 



2^ 1 „- P 2 /fe 2 



|#o(pr = (27r)--rje 

7T K Q 



(4.24) 



It is convenient to introduce a typic al diffusion coefficient 
D Q = D ko . One may use Eq. (|4.17|) to fincpl 



N / / \ 

n(r,t)=Q(t)^- t K (r/y/Dotj, (4.25) 

which decays asymptotically as n(r,t) oc exp(— r/y/Dtf) 
for r 3> y/Doi. This should be compared 



to the case where the diffusion coefficient Dq is 
momentum-independent and one finds (in 2d) n(r, t) cx 
exp (— r 2 /4:D Q t). We see, that the asymptotic profile de- 
pends crucially on the initial distribution of momenta. 
Consequently, a detailed knowledge of initial conditions 
is required for the interpretation of experiments. 



B. Diagrammatic perturbation theory for the 
nonlinear problem and the kinetic equation 

One can organize a systematic perturbation theory for 
the nonlinear problem (A ^ 0) in the limit of weak dis- 
order. This regime is characterized by the condition 
et > 1, where e is the characteristic energy determin- 
ing the diffusion coefficient. In this paper, we make use 
of the fact that in two spatial dimensions and for weak 
disorder one expects an extended regime for which the 
density evolution is diffusive, i.e. we are interested in 
nonlinear diffusion and do not consider localization ef- 
fects. We may therefore restrict ourselves to the leading 
order in the smallness paramter 1 / (er) . At this level of 
accuracy, the standard diagrammatic technique can be 
used to select diagrams for which impurity lines do not 
cross. 

In contrast to the noninteracting case discussed in the 
previous section, for which a single diffusion mode was 
sufficient for the description, the nonlinearity introduces 
an effective coupling of diffusion modes to each other. 
This coupling is not completely arbitrary, but must be 
consistent with the conservation of the total density in 
the limit of vanishing momentum. 

The relevant diagrams of perturbation theory are of 
the form depicted in Fig. |6j where the left hand side is 
associated with the initial distribution function and re- 
quires a separate consideration, see Appendix [XJ Each 
skeleton diagram, by which we mean a diagram of the 
form shown in Fig. [4j i.e., before disorder averaging, can 
be dressed by disorder in several equivalent ways, namely, 
each vertex is associated with a combinatorial factor 2. 
This is related to the fact that the interaction is chosen 
to be local in space [although we draw extended interac- 
tion lines in order to have a more convenient graphical 
representation]. This combinatorial factor is taken care 
of by choosing the decoupling in Eq. |4.34| below in two 
equivalent ways. 

We see that the expansion takes the form of a self- 
consistent Hartree theory. Due the self-consistency, the 
structure of the theory reveals itself already at the first 
order of perturbation theory in A. 

In the following we will discuss the first order perturba- 
tion theory and explain the origin of the different terms 
in the kinetic equation on this level. To this end, consider 
the diagrams in Fig. [6] The interaction line can couple 
both to the retarded and the advanced Green's function 
and due to important cancellations among these two the 
diagrams should always be grouped in pairs. The in- 
teraction line carries both momenta and frequencies and 
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FIG. 7: The two box diagrams, Br to the left and Ba to the 
right. The interaction line carries frequency uj\ and momen- 
tum qi (incoming). For the individual diagrams a constant 
term remains in the limit of vanishing external frequencies and 
momenta. This constant cancels, however, between the two 
diagrams. The cancellation is related to number conservation 
as is discussed in the main text. 



FIG. 6: On a diagrammatic level, the solution of the kinetic 
equation corresponds to the sum of all diagrams of the type 
shown in this figure. 



these can be considered as small, since they are related 
to the adjoint diffuson, or, in more physical terms, since 
the density is smooth and slowly varying in time. For 
pedagogical reasons, we will separate the discussion into 
two parts, the transfer of small momenta at vanishing 
frequency and that of small frequencies for vanishing mo- 
mentum transfer. 

We start with a finite momentum transfer. Here, the 
important point is that the diffuson to the left depends on 
the relative momentum q of the retarded and advanced 
Green's function only, but not on the sum of momenta. 
This relative momentum q is the same irrespective of 
whether the interaction line goes to the retarded or the 
advanced Green's function. 

As far as frequencies are concerned, the diffuson to 
the left of the block depends not only on the relative 
frequency, but also - via the diffusion coefficient - on 
the center of mass frequency. Therefore, it distinguishes 
between the two diagrams. 

Let us introduce the expressions for the box in the two 
cases (see also Fig. [7]) 



1 



. , (dp)G (p + -qi,E + -w 1 ) 
xG A (p_, £ _)G fl (p + , £+ ) (4.26) 



and BA(q,qi,w,a;i) = B* R (-q, -qi, -a;, -Wi). The de- 
pendence on the spectator argument e will be suppressed. 

Then the energy resolved densities for the two dia- 
grams read 

«lfl/L,e(q,w) = 

£> E (q,w) J (dqi)(dcji) n £T ^(q- q u u) - Wi) 
x0(qi,wi)Bji/i(q,w,wi) (4.27) 
where we denoted the noninteracting energy resolved 



density (compare Sec. IV A) as 

n 0jE (q, u) = 2wuf s (q, oj)V e (q, u) (4.28) 

and also used its relation to the density no(r, t) = 
J (efe)no j£ (r, i) in the linear case when introducing the 
notation 



0(r,t) = 2An o (r > *) 



(4.29) 



As will become clear in the following, i9(r, t) can be in- 
terpreted as an effective potential. 

The averaged density at order A is the sum of the two 
densities m = Uil + niR. It is 

m{q,Lj)=V e (q,u;) J {dq 1 ){dcj 1 )d{q 1 , w x ) (4.30) 

x n-o e (q - Qi, w - ui)[B R + B L ](q, q^w, w x ) 

-ujid e n 0e (q- qi,w- ui)]-[Br - B L ](q,qi,u),u>i] 

By explicit calculation one finds in the limit of small mo- 
menta and frequencies 

[Bfl + B z ](q,qi,w,wi)« -q(q-qi) (4.31) 



[Br - Bi](q,qi,w,wi) « -2i 



(4.32) 



Let us start the discussion with the case of finite mo- 
mentum transfer. Here, the combination Br + Ba enters 
the diagram and one immediately finds that the leading 
constant term cancels and the coupling is proportional 
to q(q — qi). In particular, it is proportional to the ex- 
ternal momentum q. The cancellation of the constant 
term is not accidental, but enforced by number conser- 
vation. Indeed, the limit q — >• is related to the conser- 
vation law for the total density. This can be seen best 
in the language of the kinetic equation discussed below. 
In fact, it turns out that the combination Br + Ba still 
contains a small constant term of order l/(er) 2 , which 
disappears, however, when one uses the full $ dependent 
Green's function for the self-consistent Born approxima- 
tion as is automatically the case in the kinetic equation 



approach described in Sec. IV E 
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Turning to the finite frequency transfer next, we see 
that the situation is different. Here, the constant of the 
box diagrams Br and Bl may contribute and the result is 
proportional to the difference of diffusons with different 
center of mass energies in agreement with our previous 
discussion. 

Proceedin g tow ards the kinetic equation next, one may 



multiply Eq. 4.31 by T> e (q, ui) and we present it together 



with the real space representation of the equation for no jE 



(fit - D e V 2 )n he 



-(VtfV)no £ +d t dd £ n 0e 



(d t - D e V 2 )n 0<e = 2nvf(r,t,e) 



(4.33) 



We easily recognize the first iterative solution to the ki- 
netic equation, once we use the relation between / and 
F discussed in Appendix [A] We will not follow this route 
further and formally sum up all diagrams, although this 
can be done. It has become clear that an equation is 
much more useful then any finite order in perturbation 
theory and there are more effective ways to derive the 
kinetic equation. 



C. Slow mode decomposition 

As a first step in deriving the kinetic equation we turn 
to the interaction term Si n t specified in Eq. (3.12). The 
self-consistent potential $(r, t) = 2Xn(r,t) is introduced 
in the following way. We average Si n t and obtain 



S int (4.34) 
-2A J dcdbtfat) (y>(T,t) ^(r,t)] 21 )0(r,t) 

drdt 0(r,f)#(r,t)0(r,i), 



D. Green's function 

After treating both disorder and interaction self- 
consistently as described in the previous section we ob- 
tained the action S. Due to the presence of the source 
terms describing the injection process, the fields tp and 
ip* have non-vanishing expectation values. This inconve- 
nient feature can easily be cured by shifting the fields ap- 
propriately. To this end, we introduce the Green's func- 



tion G as the average G = 



=, where the averaging 



is with respect to S = S — S s . This immediately implies 
S = f <f> G _1 0. We can define G explicitly by writing 
its inverse 

G~ 1 (r ll t 1 ,r 2 ,t 2 ) = g^ 1 {r 1 -Y 2 ,t l -t 2 ) (4.35) 

+ »7(^(ri,ti)0(r 2 ,t2))s' 

- 0(ri,ti)a(ri-r 2 )tf(ti-t 2 ). 

By denoting the averaging with the label S in this equa- 
tion, we want to remind that it should be performed with 
respect to S, not S. After introducing the shifted fields 

C( r 2,h) = 0(r 2 ,*2)-i J dr 3 G(r 2 ,t 2 ,r 3 ,0)<j) (r 3 ) 

C(ri,ti) = 0(ri,ti) + i y dr 3 o (r 3 )G(r 3 ,O,ri,ii), 



(4.36) 



we observe that S = f £ G 1 (, i-e. G — —i (CC)g- 



We 



used the fact that G 2 i = when completing the square. 
Let us also note that 



)s = (C()s + v+-F = *G + <t+-F, (4.37) 



where 



where n(r, t) — ([<j)(r, t) 
step. The averaging (. 



F(ri,r 2 ,ii,i 2 ) = 7 / dr 3 dr 4 G R (r u t x , r 3 , 0) 

Fo(r 3 ,r 4 )G A (r4,0,r 2 ,t 2 ) (4.38) 



(r,t)]i 2 ) was used in the last 
.) in both Eqs. (|4.2[) and and F (r 3 , r 4 ) = *o(r3)*o( r 4)- In particular 



(4.34) is defined self-consistently, namely wit 

S, 



respect 

Let us stress that this ap- 



tO S = S s + Sq + Sdis + Sint. 

proach includes interaction effects non-perturbatively as 
a result of self-consistency. In comparison with the clean 
case an additional factor of 2 appears in the definition 
of the self-consistent field This is not a double count- 
ing, but a result of a typical slow-mode decomposition, in 
this case in the density channel. Indeed, it will be valid 
only if $ is a slowly varying field, it means that momenta 
of the fields and <j) are close to each other. In princi- 
ple, one could also consider "anomalous" averages of the 
type (tp(r,t)tp(r,t)) and (ip* (r,t)ip* (r,t)). For systems 
for which the potential energy is not much smaller than 
the kinetic energy, such averages can in principle become 
important. In the limit we consider, namely for e 3> An, 
these terms are, however, less effective than the potential 
i? as already argued in Sec. [Tj 



n(r, t) = iG 12 (r, t, r, t) + -f(r, t, t), (4.39) 

7 



where we denoted 

/(Mi,i 2 ) =F(r,r } ^,i 2 ). 



(4.40) 



By inserting relation (4.37) into (4.35 ) we obtain an equa- 



tion for G in the form 

(id tl - ii - #(ri, tx)) G(r 1 ,t 1 ,r 2 , t 2 ) 

- J dt 3 T,(r 1 ,t 1 ,t 3 )G(r 1 ,t 3 , r 2 ,t 2 ) 

= *(n - r 2 )5(ti - h). (4.41) 

where E\ is the operator of the kinetic energy acting on 
coordinate ri. Let us comment on the different terms 
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entering the equation. The 11 and 22 components of the 
matrix G are retarded and advanced Green's functions, 
respectively, for which we use the notation G R and G . 
Disorder effects are included within the framework of the 
self-consistent Born approximation which gives rise to a 
contribution to the self-energy, 

E«K.(r,ti,t2) =7G f (r,ii,r > t 3 ), (4.42) 

The total self-energy 

£ = S dls + S s (4.43) 

additionally comprises a source term S s , which is purely 
off-diagonal and related to the initial conditions. It can 
be written as 



E s (r,ii,i 2 ) = -icr+/(r,ii,i 2 ), 



(4.44) 



wher e <j + = (c r x + i<j v )/2, and / is defined through 
Eqs. I pUO] ) and (438]). 



The equation for the Green's function (4.41) is fully 
consistent with Eq. (4.4) for the nonintcracting case. 
Here, however, G R I depend on the classical self- 
consistent potential. Besides, the dependence on the ini- 
tial conditions is explicitly included in the definition. The 
function / plays the role of the initial distribution func- 
tion in our description. The density is expressed in terms 
of the components of G as shown in Eq. (4.39). 



Thus we arrive at two equations for G and n, that are 
coupled by the self-consistency relation ■& — 2Xn. The 



first term in the Eq. (4.39) for n(r, i) accounts for diffu- 



sion for times much larger than Or, while the second 
term is a short range contribution that describes the ini- 
tial expansion up to times of the order of the scattering 
time t. It turned out to be possible to organize both the 
differential equation and the relation between the density 
n and the components of G in such a way that the in- 
formation about the initial wave function always appears 
together with G R and G A . Recall that G R and G A are 
separately averaged over disorder. 



The equation for the Green function, Eq. (4.41), still 



contains more information than is needed for calculating 
the density evolution and hence further simplifications 
can be made. In essence, we will proceed in analogy 
to the quasi-classical approximation widely used in the 
theory of nonhomogeneous superconductivity^ 1 ' ' 



E. Quasiclassical approximation 

As is well known, for the analysis of the effects of weak 
disorder, 1/t <C £ and for smooth external perturbations 
(on the scale of wave length), one may pass from the 
full quantum mechanical equations to a reduced quasi- 
classical description. In the case of superconductivity, 
this procedure leads from the Gor'kov equations to the 
Eilenberger equation in the ballistic limit and, further 
on, to the Usadel equation in the diffusive limit. Fol- 
lowing this route, we will derive an Usadel-like diffusive 



equation for a wave-packet evolving in the self-consistent 
potential which arises as a result of the nonlinearity. 
The obtained kinetic equation determines the distribu- 
tion function n(r,i,e), from which the density of the gas 
at a given moment and spatial coordinate is found as 
n(r,t) — f(ds) n(e, r, t). 

We start by introducing a mixed (Wigner) representa- 
tion for the Green's function, 

G(ri,r a ,ii,i a ) = / (dp)(de) G(r, p, t, e) 

xe ip(ri-r 2 )-ie(ti-t 2 ) ^ (445) 

where r = (ri + r 2 )/2 and t = (t\ + i 2 )/2. Considering 
first the linear case, A = 0, the frequency defines a mo- 
mentum scale p e = \/2m£, wavelength A £ = 27r/p e and 
time scale t s — e" 1 . Initially, the typical scale for e is de- 
termined by the function /(r, t, e), which in turn reflects 
the momentum distribution of the injected wave-packet, 
see Eq. (4.64) below. If the density and self-energies are 



smooth on the scale A e , the Green's function can be av- 
eraged on this scale. A necessary prerequisite is that e 
is sufficiently large. In this sense, e plays a role similar 
to the Fermi-energy in electronic systems. In the same 
spirit, the weak disorder condition, which is needed to 
formally justify the use of the self-consistent Born ap- 
proximation, can be formulated as er > 1. The averag- 
ing alluded to above can be implemented by integrating 
the Green's function in deviations from p E . In the nonlin- 
ear case A ^ 0, the frequency e in the previous argument 
should be replaced by 



e(r, t) 



0(r,i). 



(4.46) 



The quasi-classical Green's function g n can then be in- 
troduced as 

g n (r, t,s)= l -Jd£,G (r, n (p g + Vj,t,eY (4.47) 

In this equation n = p/p specifies the momentum direc- 
tion and vg — pg/m. In order to derive an equation 
for g n (r,t), one should first consider the difference of 
Eq. ( 4.41| and its conjugate equation 



G(r 1 ,t 1 ,r 2 ,t 2 ){-id t2 -i2-d(r2,t2)) (4.48) 
J dt 3 G(r 1 ,t 1 ,r 2 ,t 3 )S(r 2 ,i3,i 2 ) = S TlT2 S tlt2 . 

The result can be written as 

id t + -pv) G(r, p, t, e) - [0(r, t) • G(r, p, t, e)} 
= [E(r,f,E)'G(r,p,(, £ )]. (4.49) 
Here we introduced the »-product 

A(r,p,t,e)»B{r,p,t,e)= (4.50) 
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Due to the slowness of S and ■& a gradient expansion 
can be performed, where we keep the leading terms only. 
The quasiclassical approach in its original form does not 
involve an approximation with respect to the time argu- 
ments. Here, we make an additional smoothness assump- 
tion. Namely, we assume that the time variation of the 
density (and thereby of ■&) is sufficiently slow to justify 
the neglect of terms of the order of df-d. In addition, the 
modulus of the momentum p multiplying V is set to pi 
and the equation integrated in £, thereby obtaining an 
equation for the quasi-classical Green's function, 



id t g n (r,t,e) + -nV(p e - g a (r,t,e)) 



(4.51) 



V?9(r, t)d n g n (r, t, e) + idt'&ir, t)d e g n (r, t, e) 

Pi 

+^— [<fln(r,t,e)) n ,ff n (r,t,e)] 
= i[f(r,t,e)a + ,g n (r,t,e)} 

In this formula, ((...)) denotes angular averaging and 
d n = V n — n where V n is defined through the relation 



V p = nd p + -V„. 

P 



(4.52) 



The following relation for the disorder-part of the self- 
energy was employed 



T, dis (r,t,s) = 7 J de p v{e p ) (G(p,r,t,e)) n 



(4.53) 



-i7ri/(e)7flio(r,t,£) = -— g (r,t,e). 



The last relation serves as a definition of the scattering 
rate t^ 1 = 2irv(e)^f in our model. It was used that the 
Green's functions has a peak for e p = e, compare the 
related discussion in Sec. IIV Al 



Equation (4.51) does not fully determine the Green's 



function g a . In the quasiclassical approximation, the con- 
dition 



.9. 



^^e-^^-^^r,^) 



1. 



(4.54) 



is therefore introduced. Keeping terms that result from 
the expansion of the exponential in this formula, however, 
exceeds the accuracy of our approximation. We therefore 
use the constraint in the form 

<£(r,t,e) = l. (4.55) 

It can be seen that this constraint is consistent with 



the time evolution described by Eq. (4.51). Indeed, 



when multiplying equation (4.51 ) from the left by g n and 



adding the result to the equation that is obtained by first 



multiplying Eq. (4.51 ) by g n from the right, one obtains 



an equation for g^. The resulting equation 
d t g 2 n (r,t,e) + d t d{r 7 t)d E g 2 n {r 7 t,e) 

+v i nVgl{r,t,e) - — W(r, t)V n5n (r, t, e) 

Pi 

= -[f(r,t,e)a + ,gl(v,t,e)] (4.56) 

is solved by <7„(r,i,e) = c, where c is an arbitrary con- 
stant. This constant can be determined in the nonin- 
teracting case, where the relations g R (r,t,e) = land 
g A (r,t,e) = — 1 imply c = 1. It is usually argued ^ 16 * 39 ! 
that this constraint carries over to the interacting theory, 
and we will follow this route here. 



Equation (4.51 1, the analog of the Eilenberger equation 



in our problem, can be further simplified in the diffusive 
regime. This reduction will be discussed next. Let us 
denote by q and lo the small momenta and frequencies 
related to the space and time variation of n e . In the 
diffusive regime the inequalities rvq <C 1 and lot <C 1 
are fulfilled for typical velocities v. If we additionally 
demand rvq §/e <C 1 and lot "d/e <C 1, the main con- 
tribution comes from the zeroth angular harmonic of the 
quasiclassical Green's function. It is worth noting that 
the expansion is performed assuming that gradients and 
time derivatives of the potential are small, i.e. it is not an 
expansion in the strength of We can take into account 
the influence of higher harmonics approximately with the 
help of the ansatz 



9o + ng, 



(4.57) 



where go — (g n ) and g = d n(n'g n <) in d spatial di- 
mensions and ng is a small perturbation in the diffusive 
regime. In this limit, the constraint g\ — 1 results in the 
condition 1 = g 2 + {ng, go}, so that upon integration in 
n one obtains the relation g 2 — 1 as well as g = — <7og.9o- 
In order to derive Eq. ( 2. 1 ) , we first integrate Eq. ( 4.51 ) 



with respect to n. The result is 



id t g (r, t, e) + — V(p e ~ g(r, t, e)) 



(4.58) 



-*V0(r, t) — d - r 2 g(r, t, e) + i 9 t tf(r, t) d e g (r, t, e) 
Pi d 

= -i[/o+,0o(r,i,£)]. 



In a second step we first multiply Eq. (4.51 1 by n l before 
integrating in n and find 



id t s(r, t, e) + — V(p g g (r, t, e)) 



(4.59) 



+— W(r, t) 5o (r, t, e) + i $0(r, t) d e g{r, t, e) 

Pi 



2t>~ 



[g (r, t, e), g(r, t, e)] = -i[f<r+, g(r, t, e)] 



After multiplying this equation by g from the left and 
using the relation go[go,g] = 2g, we can formally solve 
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for g. Due to the smallness of g, not all terms need to 
be kept, and we may work with 



g(r,t,e) = g (r,t,e)V(pi g (r,t,e)) 

m 



Pi 



W(r,t, £ ) 



-hgo(r,t,e)Wg (r,t,e) 



(4.60) 



where 1^ = v^t^. We plug this expression for g into 



Eq. (4.58). In this way, we obtain the following equa- 



tion for <7o 

d t9o (r, t, e) + d t ti(r, t)d e g (r, t, e) (4.61) 
-(V + r e ~W(r, t))(D s g (r, t, e)V 5o (r, t, e)) 
= -[/(r,*,£)o-+,5o(r, *,£)], 

where = VgTg/d. We used the relation ie/pg x (2 — 
d)/d — TgDg, where we defined the quantity 



-c? £ lnv(e). 



(4.62) 



r vanishes in two spatial dimensions, since the density of 
states is constant. In three dimensions, however, T = 
is finite. As mentioned before, the above equation should 
be supplemented with the matrix constraint <7g(r, t,e) = 
1. 

Before making a specific ansatz for the solution, let 
us focus on the function / that specifies the injection of 
the wave-packet and initial evolution up to times of the 
order of the scattering time r, see Eq. (4.40). If F(p, r) 



is sufficiently smooth in the sense that for typical v = 
p/m and q controlled by F(p, q) = ^o(p + q/2)*g(p - 
q/2) the inequality rvq <C 1 holds, we can approximately 
replace 

2-kv s f{v,t,e) (4.63) 
« 5{t) y"(dp)F(p,r)27n5( £p + tf(r,0)-£) 
= 8{t) F(£-tf(r,0),r). 
Next we introduce the following ansatz for g : 



—^-hir, t, e) 



(4.64) 



which solves the equation provided n fulfills the kinetic 
equation 

d t h(r, t, e) - V(D g Vrh(r, t, e)) 
+d t $(r,t)d £ n(r,t,e) = 6{t)2TT^(e)F(e, r). (4.65) 

Here, we used the notation Vr = V — Vz?(r, t)Tg. This 
concludes our derivation of the kinetic equation from 
Eqs.(4.41). The diagrammatic interpretation of the dif- 



ferent terms appearing in this equation was provided in 
Sec. IIVBI for the two-dimensional case. The main new 
ingredient for d ^ 2 is the non-constant density of states. 
Within our model, it results in a frequency-dependent 



scattering time [compare Eqs. (4.62)]. Since the density 
of states enters with argument e = e — $(r, t), the disor- 
der part of the self-energy E^ is explicitly depends on 
In a diagrammatic language, it means that a generaliza- 
tion of the box diagrams B is required for a non-constant 
density of states in order to accommodate this change, 
see Fig. [8] This modification was first noticed in Ref. ITS! 





FIG. 8: For dimension d 7^ 2 the density of states is not con- 
stant and the ^-dependence of the Green's function entering 
the SCBA becomes important. In this case the box diagrams 
Br and Ba should be generalized as displayed. 

One may write the distribution function as a function 
of the kinetic energy instead of the total one. n(r, s,t) = 
n(r, e + 'i?(r, t), t), see Eq. (2.5). From a technical point of 



view, this transformation amounts to a gauge transfor- 
mation. We could have utilized this transformation al- 
ready at the beginning of our derivation by working with 
the so-called gauge- invariant Green's function Q, which 
can be introduced as 

G( ri , r 2 , i l5 i 2 ) - J (dp)(de) Q(r, p, t, e) 

e »p(ri-r 2 )-i[e— #(r,t)](ti-t 2 ) /4_gg\ 

A derivation based on Q instead of G, but otherwise fol- 
lowing the same lines as described in this section, leads 
directly to Eq. (pl^ instead of Eq. J2~T). 



V. COLLISIONS INDUCED BY THE 
NONLINEARITY 

The purpose of the paper is to present the technical 
aspects of the derivation of the kinetic equation describ- 
ing the pulse propagation in a disordered and nonlin- 
ear medium. The obtained equation (2.5) describes the 



diffuse propagation (as a result of collisions with elastic 
defects) in the self-consistent potential, but so far fully 
ignores collisions induced by the nonlinearity. With re- 
spect to the nonlinearity, this equation describes the col- 
lisionlcss regime. 

We now wish to discuss the role of collisions. To this 
end let us first recall the general spirit of the deriva- 
tion presented for the collisionsless regime. As is typical 
for disordered systems, the physics at long time scales 
and long distances is dominated by diffusion modes. The 
nonlinearity leads to an interaction of these modes. To 
treat this effect, we singled out pairs of fields 4> and 
d> with a small momentum difference in the interaction 
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term. Afterwards, the effect of interaction of the diffu- 
sion modes was considered in a self-consistent way by 
introducing the smooth classical potential •& as described 
in Sec. |1V C| This procedure may be viewed as a mean 
field approximation, ft should be noted, however, that 
in this procedure only a small (albeit important) subset 
of all possible scattering processes was singled out and 
treated non-perturbatively as a result of self-consistency. 
It is important that the potential is proportional to the 
density and therefore smooth and slowly varying. This 
is the reason why the self-consistent part of the problem 
of the propagation of the diffusion modes may be treated 
within the quasiclassical formalism. 

To incorporate collisions, one has to go beyond the 
scheme discussed above. Collisions induced by a non- 
linear interaction in classical wave systems are routinely 
studied in nonlinear physics (see e.g., V. E. Zakharov, 
V. S. L'vov, and G. Falkovich "Kolmogorov Spectra of 
Turbulence"). There it works as follows. In the equa- 
tion of motion for the occupation numbers n p (r, t), one 
obtains nonlinear terms, which are considered using the 
random-phase approximation. At second order in the 
coupling constant A this procedure yields a collision inte- 
gral, which in the case of the four-wave interaction (like 
in the NLSE or GPE) is proportional to the third power 
in the occupation numbers. 

Here we will show how the derivation of the collision 
integral in the kinetic equation can be obtained in the 
framework of the field-theoretical approach we use. In 
order to account for collisions, we have to go beyond 
the mean-field description employed in the collisionless 
regime, namely, we need to include fluctuations. We will 
derive these at the second order with respect to A, the 
lowest order at which collisions appear in the theory. We 
therefore need to calculate second-order corrections to 
the self energy. When doing so, we will assume that the 
diffusive propagation in the field of the self-consistent 
potential created by the nonlinearity is already known 
according to the analysis presented in the previous sec- 
tion. 

In the calculation of the collision integral, we will 
use the G reen functions G(ri,t 1 ;r 2 ,t 2 ) as defined in 
Eq. (4.35). When doing so, we neglect terms contain- 
ing F, cf. Eq. (4.37), because F decays on a scale of 
the mean-free path in the disordered medium. The off- 
diagonal component G\ 2 , in turn, describes the long- 
range nonlinear diffusion of a partial wave until the mo- 
ment of collision with another partial wave at t t\ w t 2 . 
The component G\ 2 resembles the Keldysh component in 
the regular technique. It is non-vanishing due to the in- 
jection process at t = 0, which is encoded in the source 
term S s of the action. We will therefore denote Gyx as 
G s . 

Since Si n t originates from the NLSE/GPE, the theory 
used in this paper contains only classical vertices, namely 
those that couple one of the quantum components of the 
doublets or with three classical ones, compare Fig. [3] 
As a consequence of this fact, there are only three con- 




FIG. 9: Corrections to £ B according to Eq. 



5.2 




FIG. 10: Corrections to £ according to Eq. 5.1 



tributions (diagrams) to be calculated for self-energies: 
one for S s and the other two for each of the diagonal 
components, e.g., for T, R . As a result one obtains 



S 5 (xi,x 2 ) = -2X 2 G S \x 2 , Xl )G b ; ( Xl ,x 2 )G b \x u x 2 ), 



and 



(5.1) 



E R (x u X2) = -2\ 2 [G A {x 2 ,x 1 )G s {x l ,x 2 )G s {x u x 2 ) 
+2G s (x 2 ,x 1 )G R (x 1 ,x 2 )G s (x 1 ,x 2 )}. (5.2) 

These two quantities determine the collision integral 
in the kinetic equation. In standard kinetic theory, the 
kinetic equation is formulated in terms of the mass-shell 
distribution function n p (r, t). It can be introduced as 
follows. First, one parametrizes G s = G R • n — n • G A , 
where n = h(r, p, t, e) and the "-product has been defined 



in Eq. (|4.50|). Then one defines the on-shell distribution 
n(r,p,t,s = e p + 'd + $t[E R )) 



function as 
n p (r,t) 



(5.3) 



This definition is motivated by the observation that as 
long as G s is a smooth function of coordinates and times, 
the largest contribution to G s comes from the product 
of Wigner transforms 

G s (r,p,t,e) » n(r,p,t,e)(G R -G A )(r,p,t,s) (5.4) 



-2win(r, p, t, e)5(e 



i9 - M(T, R )) 



re -27rin p (r,t)5{e-e p -$-Vi(Y, R )). 

Here, 5 is a broadened 5-function, which is sharply 
peaked compared to the scale of variation of n and can 
therefore be replaced by a regular delta function. Since 
the distribution function n always appears in combina- 
tion with the delta function, it is useful to work with the 



mass-shell distribution function n p defined in Eq. (5.3) 
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In a general context, could be an external poten- 
tial. For our application, an external potential is not 
present, but we treat a part of 3?E fl separately, namely 
the self-cons isten t potential i9(r, t) = 2An(r, t). There- 
fore, in Eq. (5.4 1 should be understood as the real 
part of the E K as given in Eq. |9j). As will be discussed 
further below, in the regime of applicability of our ap- 
proach may be considered to be small compared to 

and we will not mention it further. 

In an approximation consistent with this reasoning, the 
collision integral can be written as 

F M (r,p,t,s) (5.5) 
= iS s (r, p, t, s) + 2n(r, p, t, e)$SH H (r, p, t, e). 



Entering with Eq. (5.4) into the expressions (5.1) and 



( 5.2 ) , and replacing the Wigner transform of the products 



of Green's functions on the RHS of both equations by the 
product of Wigner transforms, one finds from Eq. (5.5) 



the collision integral for the mass-shell distribution func- 
tion 



7 coll (r,p,i) = I(r,p,t,e = e p + #) 
dp 2 dp 3 ap4 



(5.6) 



2,/ 



(27T) 

S(p + p 2 - p 3 - Pi)S(e p + e P2 - e P3 - £ P4 ) 
x { [n p + n P2 } n P3 n P4 - n p n P2 [n P3 + n P4 ] } , 

where we suppressed (for the sake of brevity) the space 
and time arguments x — (r, t) in the distribution func- 
tions. The normalization is n(r,t) — J (dp)n p (r, t). One 
may check that the obtained collision term coincides with 
the C22-term in the kinetic theory of the Bose gasp^the 
Bose- factors n p + 1 are replaced by n p . Most important 
is that the obtained collision integral is universal, i.e., it 
does not depend on the system microscopy and holds for 
both NLS and GP equations. 

In this paper, we work with the ^-integrated quasi- 
classical Green's function and, correspondingly, with a 
frequency-dependent distribution function rather than 
with a distribution function depending on the quasipar- 
ticle energy e(p). To make a con nection, we note that 
in the equation for G s , Eq. (5.4), where the <5-function 



was used to fix the frequency argument e of n, we may 
alternatively fix e p and thereby the modulus of p 



G s (r, p, t, e) « —2irih(T,pgn,t,e)6(i— e f 



(5.7) 



where e is defined in Eq. (4.46). With the help of this 
representation one finds 

C n (r, t, e) = I cM (r, p E n, t, e + 0(r, t)) (5.8) 

= 47rA 2 (27r) d J dn.2dn.3dn4 J ds2desds4 

y.v{e% s )v{e^)v{ei)^{e + e 2 - e 3 - e 4 ) 

X S(p e + p e2 - p E3 - p E4 ) ( [n' + nl 



where p £i = p ei n and we suppressed the space and 
time arguments in the distribution functions n' n e (r, t) = 
n(r,p e n, t,e + i?(r, t)). Only positive values of Si are in- 
cluded in the integration. 

In the diffusive limit, which we concentrate on in this 
paper, only the isotropic part of n' a is important and n' a 
may be replaced by its angular average n' = J dn n' n . 



In a similar way, the knowledge of T 



coll 



Jdn I 



coll 



sufficient. In accordance with Ref. I13| the normalization 
of the distribution function n(r, t, e) in Sec. |iVE"| has been 
chosen such that n(e) = 2 / nv{e)n' '(e). The resulting full 
kinetic equation including the collision integral is written 
in Sec.|nJ Eq. fl2T3 ). 

In the course of derivation of the kinetic equation we 
omitted the renormalization induced by the real part of 
Yi R . This kind of renormalization is standard for any 
many-body problem. The corrections induced by the real 
part, for example 9 £ !R£ , are of the order of (An/e) 2 , 
and are therefore smaller than the leading terms which 
are kept in the kinetic equation. 



VI. CONCLUSION 

In this work, we discussed the propagation of a wave- 
packet in a disordered and nonlinear medium, for which 
the dynamics is governed by the NLSE/GPE. Possible 
applications include the propagation of a light beam in a 
nonlinear optical medium and the expansion of a cloud 
of Bose atoms released from a trap. For dcfinitcness, we 
use the term "particles" irrespective of the system. 

We considered the case when the potential (interac- 
tion) energy induced by the nonlinearity is considerably 
smaller than the typical kinetic energy This allowed us 
to use the picture of a gas of particles moving in a self- 
consistent potential rather than that of a hydrodynamic 
flow. Diffusion occurs as a result of elastic scattering from 
a random potential. We studied a regime for which par- 
ticles scatter on impurities many times before colliding 
with other particles. 

Another important consequence of the smallness of 
the nonlinearity is the possibility to neglect off-diagonal 
terms in the Bogoliubov transformation. In the case of 
a Bose condensate released from a trap, our considera- 
tion corresponds to a stage of evolution when the ini- 
tial hydrodynamic fk W 40 * 41 ! of the Bose atoms already 
passed by and particles diffuse with a typical kinetic en- 
ergy of the order of the (initial) chemical potential and 
the wavelength \ typ comparable with the healing length £ 
of the trapped condensate. We assume that Xt yp is much 
shorter than the mean free path. Since we use the GPE, 
which arises as the classical equation of motion in the 
theory of the interacting Bose gas, it is assumed that the 
occupation numbers n p with p ~ 2-7T /Xtyp remain large 
on the discussed stage of the expansion. 

Compared to the case of disordered electrons with 
electron-electron interactional virtual processes involv- 
ing diffusion modes need not be considered for the dis- 
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cussed problem. As we have already mentioned, such 
processes give corrections that are small in the parame- 
ter 1/Et -C 1, but unlike for electrons at low temperature, 
they are not accompanied with non-analytic corrections, 
which make them important in the case of the degenerate 
electron gas. 

In the case of two-dimensional particles, d = 2, the mo- 
tion in the plane is not constraint, while the third dimen- 
sion cither represents the effective time-like direction in 
the case of optics experiments or is blocked by the quan- 
tization induced by a potential that restricts the motion 
in the transverse direction. It is important to distinguish 
the original dimension of the single particle states in the 
NLSE/GPE, denoted with d in this paper, from the ef- 
fective dimension of the diffusive collective modes, which 
may be different. Namely, the derived kinetic equation 
can be solved in different geometries. 

To illustrate the role of the effective dimensionality, 
let us consider the example of a stripe made out of 
2d-particles. Then, the diffusion will be described by 
a one-dimensional solution while particles can be two- 
dimensional if the quantization with respect to the width 
of the stripe can be neglected. As in the case of 2d- 
particles diffusing in a plane, the exact solution for the 
time-dependence of the mean squared radius still holds. 
We expect that the existence of this simple analytical 
result can be useful for numerics or suitably designed ex- 
periments. 

Besides technical details of the derivation of the 
two self-consistent equations describing the collisionless 
regime, the present paper contains a discussion of inter- 
particle collisions. The collision integral has been ob- 
tained from the same field-theoretical approach that was 
used as a starting point for the derivation of the kinetic 
equation in the collisionless regime and the procedure 
was straightforward. It is important to stress that the 
collision integral is the same for both optics (NLSE) and 
cold atoms (GPE), i.e., independent of the microscopic 
origin. Since the inter-particle collis ions are elastic and 
local, it does not alter the relation (2.12) between (r 2 ) 
and t in the case of a constant density of states. The 
result remains valid as long as the kinetic equation is ap- 
plicable, despite the fact that the collisions change the 
dynamics of propagation. 

The change of dynamics caused by collisions is qual- 
itatively different from the one introduced by the self- 
consistent potential. Indeed, the smooth self-consistent 
potential is responsible for a gradual change of the kinetic 
energy during the expansion. In contrast, the energies of 
incoming and outgoing particles participating in a col- 
lision process may differ considerably. In particular, in 
three dimensions the collision-induced redistribution of 
energies may lead to a population of localized particles 
with energies e < 1/r. This mechanism bears a certain 
similarity with the seeding of a macroscopic occupation 
of low-energy states in a trapped Bose gas; this step is 
crucial for the formation of a Bose-condensate starting 
from a confined Bose gas^HUl. I n the case of the expand- 



ing disordered Bose gas in 3d, collisions seed a population 
of particles that are likely to localize. An estimate for the 
rate of generation of localized particles may be obtained 
from the in-scattering term of the collision integral upon 
integration over the interval < e < 1/r, namely 

, ,,, 1 n'(e ~ 2e) n(t) 

dn loc /dt « dn/dt^ 1/T ~ — -i— ^^ -(6.1) 

Here, we assumed that both colliding particles have the 
energy ~ e. Although the discussed effect is important 
under static conditions, one may show that for the sit- 
uation we study the seeding of localized states for an 
expanding cloud is negligible because of the fast drop of 
n(t). 

The situation in two dimensions is different in that in 
the absence of interactions all states are localized on the 
scale of the localization length, ^ oc (e). For a state with 
the energy e, the process of localization starts to develop 
at a time of the order of lf oc (e)/D(e). In the presence of 
interactions, however, this picture changes. First of all, 
both the time-dependent potential and the interparticle 
collisions lead to dephasing, which weakens localization 
effects. One may show that in 2d if the number of par- 
ticles is large enough, the expanding cloud will pass k oc 
without being stopped. We will discuss this situation in 
more detail elsewhere. 
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Appendix A: Distribution function 

The aim of this appendix is to show how the approxi- 
mation (4.22) for the distribution function / is obtained. 



For the sake of completeness, we consider the general- 
ization to the case with interaction. Starting point is 
the definition of / in Eq. (4.18). We introduce times 
t = (ti + t2)/2 and At = ti — t% as well as coordinates 
r = (ri + r 2 )/2 and p — T\ — r 2 and write 

G(r,p,i,e) = y"d(Ai)dpG(r 1 ,r 2 ,i 1 ,i 2 )e- 4p ' ,+l£A '. 

(Al) 

Note that the Green's function G depends on coordinates 
r, t only via the self-consistent potential $, since other- 
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wise the disorder averaged system would be translation- 
ally invariant in time and coordinates. If variations of i? 
in time in space are slow in comparison to other relevant 
scales in the system, we may write the leading term in a 
gradient expansion for /(r, t, e) (defined in Eq. (4.21 )) as 



/(r,M) « 7 J(dp)(dq)(diu)G R (v,p +1 t,e + ) 

xF(p, q)e iqr G A (r, p_, i, e_) e~ luJt .(A2) 

Both Green's functions decay on typical time scales of the 
order of the mean free path. If / is convoluted with a 
function that is smooth of this time scale, we may there- 
fore use f(r,t,e) ~ 8(t) f_ dtf(r, e, t) instead. Next, it 
is assumed that F controls momenta q so that essential 
q = |q| are small as qlt yp <C 1- Then, 



J(r,t,e)*6{t) /(dp)(dq) 



7 F(p,q)e^ r 



•0(r,O)] s 



(2r f ) 2 
(A3) 



The Lorentzian is peaked around e ~ e p + i?(r, 0) and 
has a width of the order of t. 

If the distribution function is used to average a quan- 
tity that depends smoothly on e as is the case for our 
problem (where e determines the diffusion coefficient), 
the Lorentzian acts essentially like a smeared S func- 
tion and we may use f(r,t,e) = Tg(2n)5(e — e p — 



0(r,O))/(<fe)/(r,t, E ) 
appendix, 



This leads us to the result of this 



2nv e f(r,t,e)= (A4) 
5(0 / (dp) F(p, r) (2tt)J (e p + 0(r, 0) - e) 



We used the relation jt§ = l/(27ri/g). If the smoothness 
assumptions outlined in this appendix are met, the repre- 
sentation of the distribution function / in the form given 

by / is justified. In the noninteracting limit, this leads 
us to relation (14. 221). 
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